Cracks and Crazes: On calculating the macroscopic fracture energy of glassy polymers 

from molecular simulations 
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We combine molecular dynamics simulations of deformation at the submicron scale with a simple 
continuum fracture mechanics model for the onset of crack propagation to calculate the macroscopic 
fracture energy of amorphous glassy polymers. Key ingredients in this multiscale approach are 
the elastic properties of polymer crazes and the stress at which craze fibrils fail through chain 
pullout or scission. Our results are in quantitative agreement with dimensionless ratios that describe 
experimental polymers and their variation with temperature, polymer length and polymer rigidity. 

PACS numbers: 82.35.Lr, 83.60.Uv, 62.20.Mk 
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Understanding the molecular origins of macroscopic 
mechanical properties such as the fracture energy G c is a 
fundamental scientific challenge. In tough materials, the 
work G c required to propagate a crack through a unit 
area is orders of magnitude higher than the lower bound 
provided by the equilibrium interfacial free energy G eq 
of the crack surfaces. Efforts to calculate this large in- 
crease in fracture energy have been frustrated, because 
phenomena on many length scales must be treated simul- 
taneously [1] . In both amorphous and crystalline materi- 
als, the fracture energy depends on processes that range 
from breaking of atomic bonds to formation of defect 
structures on micron and larger scales. 

In this Letter, we present a multiscale approach that 
allows us to calculate the plane strain fracture energy of 
an important class of unfilled amorphous polymer glasses. 
Experiments [4-6] show that under tensile loading the 
fracture energy of materials such as polystyrene (PS) and 
polymethylmethacrylate (PMMA) is mainly due to the 
formation of an intriguing craze structure in a "process 
zone" around the crack tip (Fig. 1). This results in a 
large increase in fracture energy, G c /G eq ~ 10 3 — 10 4 , 
that is essential to the use of these materials as adhe- 
sives, packaging materials and windows [2-4]. 

In the craze, ~ 0.5nm diameter polymers are bundled 
into an intricate network of ~ lOnm diameter polymers 
that extends ~ 10[im on either side of the ~ mm crack 
and ~ 100/zm ahead of the crack tip. Molecular level sim- 
ulations of regions with linear dimensions of mm or even 
jim are not feasible. They would also be inefficient, since 
most regions near the crack are homogeneous enough to 
be treated with continuum mechanics [5,7]. Here, we 
combine the two approaches using molecular simulations 
of representative volume elements (see Fig. I) to provide 
information about craze formation, deformation and fail- 
ure that is needed to construct a continuum fracture me- 
chanics model. 

One advantage of studying polymeric systems is that 
many dimensionless ratios are independent of the specific 
chemistry of the molecules. For this reason we consider 



a bead-spring model that has been shown to provide a 
realistic description of polymer behavior [8-11]. Each 
linear polymer contains N beads of mass m. Van der 
Waals interactions between beads separated by a distance 
r are modeled with a truncated Lennard-Jones potential: 
V LJ (r) = Au [(a/r) 12 - (a/rf - (a/r c ) 12 + (a/r c ) 6 ] for 
r < r c = 1.5 a, where uq and a are characteristic en- 
ergy and length scales. A simple analytic potential, 



Vb r (r) 



-ki(r — r c ) s (r — R\), is used for the cova- 



lent bonds between adjacent beads along each chain. 
The constants k\ and R± are adjusted to fix the equi- 
librium bond length [8], 0.96a, and the ratio of the 
forces at which covalent and van der Waals bonds break. 
We find that this ratio is the only important param- 
eter in the covalent potential and set it to 100 based 
on data for real polymers [11,12]. The polymer rigid- 
ity and entanglement length iV e are varied by intro- 
ducing local bond-bending forces [10] with a potential 
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bone, fi denotes the position of the zth bead along the 
chain, and b characterizes the stiffness. Two limiting 
cases of fully flexible (N e ~ 60 — 70 beads, b = 0u ) and 
semiflexible (N e ~ 30 beads, b — l.Stto) chains are con- 
sidered here. The chain length is varied from N = 64 
beads to 1024 beads. 

We first show that our model captures the essential 
experimental features of craze formation (Fig. 1 A) . The 
simulation cell has periodic boundary conditions and is 
initially a cube of size L. The length along one direction 
L3 is increased at constant rate, while the other dimen- 
sions of the cell are held fixed. Fig. 2(a) shows typical 
results for the stress 03 along the stretching direction as 
a function of the elongation L3/L. In all cases, there is 
an initial peak at small strains, where the material yields 
by cavitation [13,14]. As in experiment [6], this peak is 
followed by a long plateau at a constant stress S. During 
this plateau, deformation is localized in a narrow "ac- 
tive zone" at the boundary of the growing craze network 
(Fig. 1A). 
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FIG. 1. Some length scales involved in the fracture of glassy polymers. A ~mm long crack penetrates the material. A craze of 
typical width d — 10 — 50/xm and length I ~ 10<f forms in the "process zone" in front of the crack tip (top left). Craze formation 
(A), deformation (B) and failure (C) are studied with molecular dynamics simulations of ~100 nm-sized volume elements. A: 
craze fibrils emerge from the dense polymer glass (top) within a narrow "active zone", B: fully developed craze structure, C: 
failure of the craze immediately in front of the crack tip. Individual beads of the polymer chains are shown as small ellipsoids. 
Chains carrying the largest tension are colored blue and broken bonds in C are colored red. The width of each image is 128 
bead diameters. The fibril spacing Do is ~10 bead diameters for this figure, but depends on temperature, chain rigidity and 
other parameters. 



S represents the stress needed to draw fibrils out of the 
dense regions adjacent to the craze [9]. This steady state 
drawing process increases the volume occupied by the 
polymer by a constant factor called the extension ratio 
A. When L3/L reaches A, the entire system has evolved 
into a craze, and the stress (Fig. 2(a)) begins to rise. 

It is evident from Fig. 2(a) that A is strongly depen- 
dent on chain rigidity and thus the entanglement length 
N e . We find that A decreases from about 6.1 for flex- 
ible chains to 3.6 for semiflexible chains. As in exper- 
iments, these values are quantitatively consistent with 
a simple model that assumes entanglements act like per- 
manent chemical crosslinks [6]. During crazing, segments 
between entanglements are expanded from their equilib- 
rium random-walk configurations to nearly straight lines 
(Fig. IB). 

We consider the common case where the dominant con- 
tribution to the fracture energy is the work needed to 
craze material in the process zone ahead of the crack tip 
[4] . As the crack advances, each region is expanded 



at the constant plateau stress S. In steady state, ad- 
vancing the crack over an area A has the net effect of 
expanding a region of this area at constant stress from 
its initial width to the final width d at which the craze 
cracks. Thus G c — S(d — d/X) [15]. After normalizing 
by the lower bound for the fracture energy provided by 
the interfacial free energy change G eq = 27, the fracture 
energy can be written in dimensionless form as 
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(1 - 1/A). 
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Eq. (1) shows that G c is primarily limited by the craze 
width d. Unlike the other quantities in Eq. (1), d cannot 
be obtained directly from MD simulations. However, a 
minimal continuum model proposed by Brown [5] allows 
us to calculate d. 

Brown pointed out that although there is a constant 
plateau stress S on the craze boundary, a stress con- 
centration occurs near the crack tip (see Fig. 1). He 
formulated a fracture criterion by equating the stress at 
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the crack tip to the maximum stress S max that the craze 
fibrils can withstand. The stress at a distance r from 
a crack tip in a continuous elastic medium diverges as 
r -i/2 This divergence is cut off at the characteristic fib- 
ril spacing Dq, below which the material can no longer 
be treated as a homogeneous elastic medium. Since the 
stress varies from SVax to S as r varies from Do to d/2, 

1 /2 

Smax oc S (d/ Dq) . Solving the linear fracture mechan- 
ics problem within the craze yields [16] 
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where S mayi denotes the maximum stress that the craze 
fibrils can withstand, and the prefactor k depends on 
the anisotropic elastic constants of the craze network 
(Eq. (3)). Neither the elastic properties of the craze 
network nor S max are easily obtained from experiments. 
However, we can calculate both from MD simulations of 
regions B and C in Fig. 1 and thereby provide the key 
ingredients for calculating the fracture energy of glassy 
polymers from Eqs. (1) and (2). 

Elastic constants were calculated by applying small 
(< 0.5%) step strains to fully developed crazes (Fig. 
IB) at two different elongations L 3 /L and measuring the 
change in stress. Table I shows key ratios for flexible and 
scmiflexiblc chains at two representative temperatures 
T = 0.1 u /ks and T — O.Oluo/fcs. Both are well below 
the glass transition temperature T g « 0.35wo/fc_B °f tnc 
bead-spring model. As can be expected from the highly 
oriented structure of the craze network (see Fig. 1), C33 
is always much bigger than the other elastic constants, 
which are all of the same order. The prefactor k in Eq. (2) 
is given by [16] 



(1-C2) + (C33/2C44)(1-Ci) 

2(1 -C1) 2 



(3) 



where C\ = 2C2C13/C33 and C 2 = 013/(011 + 012)- Insert- 
ing the elastic constants, we obtain values for k between 
2.0 and 2.8 for flexible chains and between 1.1 and 1.7 for 
semiflexible chains. The crazes with higher elongations 
always have a lower value of k. 

A simple approximate expression k « \J C33 /4c44 can 
be obtained by noting that C33 3> C13, and thus C\ ~ 
and C2 ~ 1. Table I also shows that this is an accurate 
approximation for all practical purposes. This simple 
expression shows clearly that the ability of crazes to re- 
sist shear (C44 > 0) limits their fracture energy. As first 
pointed out by Brown [5], the absence of lateral stress 
transfer would lead to n — > 00 and thus to an infinite G c . 

To determine the stress S maK at which fibrils break, we 
continue straining the fully developed craze until it fails 
(Fig. 1C). Although all chains that are long enough to 
form stable crazes (N/N e > 2 [9]) show the same plateau 
stress and extension ratio, their crazes exhibit very dif- 
ferent behavior for L3/L > A (Fig. 2(a)). Short chains of 



length N/N e — 2 easily pull free from the topological con- 
straints imposed by entanglements, and the stress drops 
monotonically. As N increases, the force needed to pull 
chains free from entanglements along a failure plane rises, 
and there is a corresponding increase in S max . The failure 
mechanism changes when the force needed to disentangle 
the chains reaches the breaking force for covalent bonds. 
At this point the forces along the chains and S max both 
saturate due to chain scission. 
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FIG. 2. (a) Normalized stress 03/ S during craze growth 
at T = 0.3uo/fcs for flexible (solid lines) and semiflexible 
(dotted lines) chains of the indicated length in systems com- 
posed of 32768 beads. Here we take N e = 32 for the semi- 
flexible and iV e = 64 for the flexible chains. The value of 
Smax/S 1 is obtained from the maximum height of the curves. 
Panel (b) summarizes Sma^/S for the flexible (o) as well as 
semiflexible (A) chains at T = 0.3iio/fcs (open symbols) and 
T — 0.1uo/Jcb (filled symbols) as a function of chain length. 

Fig. 2(b) summarizes our results for S max /S as a func- 
tion of chain length, temperature and flexibility As 
N/N e rises above 2, 5 m ax/<5 rises rapidly and then sat- 
urates due to the change in failure mechanism from 
chain pullout to chain scission. Saturation occurs for 
N/N e between 8 and 16. The limiting value of S max /S 
lies between 3.4-3.8 for T < 0Au /k B and 5.0-5.3 for 
T = Q-Zuo/ks and increases slightly with chain rigidity. 
In general, more chain scission is observed for chains with 
higher rigidity and at lower temperatures. 

We are now in a position to evaluate Eq. (2) and 
compare our results to experimental values. For flex- 
ible chains, k = 2.0 - 2.8 and S max /S = 3.4 - 5.0, 
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yielding d/Do between 290 - 890. Values for semiflexible 
chains give d/D between 200 - 600. Typical craze widths 
d observed in experiments range between 3 — 20 fim, 
whereas characteristic fibril spacings Dq have values be- 
tween 20 — 30 nm. Thus the range of experimental values 
for d/D = 100 — 1000 overlaps well with our results. 

In addition to values quoted above, calculating the 
fracture energy from Eq. (1) requires values for the 
plateau stress S, mean fibril spacing D , and surface 
tension 7. Typical values from our simulations are S = 
0.5- 1.4u /a 3 , D = 10- 14a, and 7 = 0.6 - 1.0 u /a 2 . 
With these values, we arrive at our final result G c /G cq = 
1300-4300 (flexible polymers) and G c /G cq = 1200-3500 
(semiflexible polymers). Within this range, G c /G eq tends 
to drop with increasing chain rigidity and decreasing tem- 
perature T. This trend is also found in real adhesive 
joints, where the fracture energy generally decreases with 
decreasing temperature [2]. 

Our simulations agree with experimental observations 
in greater detail. Sha et. al. [17] have compiled values of 
G c for PS and PMMA as a function of polymer molecu- 
lar weight M w . Neither polymer shows a large fracture 
energy when M w is less than twice M e . As in our sim- 
ulations (Fig. 2(b)), the fracture energy rises rapidly as 
M w /M e rises above 2 and then saturates around 10 M e . 
The limiting values of G c /G eq at large M w /M e , 2500 
(PMMA) and 5000 (PS), are comparable to our predicted 
values. 

In conclusion, we have demonstrated that supplement- 
ing a simple continuum model with constitutive relations 
from molecular simulations can provide quantitative pre- 
dictions for key material parameters such as the fracture 
energy. This approach can be further developed by us- 
ing chemically detailed interaction potentials for specific 
polymers in the molecular simulations. For example, one 
might expect that a realistic treatment of side groups 
could increase the friction between polymers during craze 
formation, and increase the likelihood of chain scission. 
A more detailed finite element model for crack propaga- 
tion could also be used, as in recent work [18,19] that 
assumed simple constitutive relations for craze widening 
and breakdown. Although one might hope to include all 
length scales simultaneously in a hybrid calculation, this 
is complicated by the rapid increase in the relevant time 
scale with increasing length scale [7] . Finally, it should be 
noted that macroscopic cracks can contain an ensemble 
of crazes with characteristic sizes and spacings and that 
other processes such as shear banding can contribute to 
the fracture energy, depending on loading conditions and 
materials. These issues should provide fruitful topics for 
future work. 
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TABLE I. Elastic properties of model crazes composed of 
262 144 beads with flexible (fl.) and semiflexible (sfi.) chains 
at two elongations L3/L. Uncertainties in dj are about 10%. 
Here N = 256, but the results do not depend on N for 
N > 2N e . They are also insensitive to the covalent bond 
potential, since strain is accommodated by the weaker van 
der Waals bonds. 





T[u /k B ] 


L 3 /L 


C11/C33 


C44 / C33 


K 


\J C33/4C44 


fl. 


0.01 


5.5 


0.026 


0.038 


2.8 


2.6 




0.01 


7.9 


0.016 


0.065 


2.0 


2.0 




0.1 


5.8 


0.030 


0.041 


2.7 


2.5 




0.1 


7.9 


0.015 


0.054 


2.2 


2.1 


sfl. 


0.01 


3.4 


0.12 


0.10 


1.7 


1.6 




0.01 


4.8 


0.051 


0.10 


1.6 


1.6 




0.1 


3.4 


0.087 


0.086 


1.9 


1.7 




0.1 


4.8 


0.026 


0.15 


1.4 


1.3 
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